const std = @import("std");
const math = std.math;
const lbmath = @import("./math/lbmath.zig");
const Allocator = std.mem.Allocator;
const dprint = std.debug.print;

inline fn type_of_matrix(comptime T: type) u32 {
    const info = @typeInfo(T);
    switch (info) {
        .array => |a| {
            switch (@typeInfo(a.child)) {
                .array => return 1,
                else => return 0,
            }
        },
        .pointer => |p| {
            switch (p.size) {
                .one => {
                    switch (@typeInfo(p.child)) {
                        .array => |a| {
                            switch (@typeInfo(a.child)) {
                                .array => return 2,
                                else => return 0,
                            }
                        },
                        else => return 0,
                    }
                },
                .slice => {
                    switch (@typeInfo(p.child)) {
                        .pointer => |pp| {
                            switch (pp.size) {
                                .slice => return 3,
                                else => return 0,
                            }
                        },
                        else => return 0,
                    }
                },
                else => return 0,
            }
        },
        else => return 0,
    }
    // return re;
}
inline fn type_is_vector(comptime T: type) bool {
    const info1 = @typeInfo(T);
    switch (info1) {
        .pointer => |pp| {
            switch (pp.size) {
                .slice => return (@typeInfo(pp.child) != .pointer and @typeInfo(pp.child) != .array),
                .one => return (@typeInfo(pp.child) == .array and @typeInfo(@typeInfo(pp.child).array.child) != .array),
                else => return false,
            }
        },
        .array => |a| return @typeInfo(a.child) != .pointer and @typeInfo(a.child) != .array,
        else => return false,
    }
    // return info1 == .pointer and info1.pointer.size == .slice and @typeInfo(info1.pointer.child) != .pointer and @typeInfo(info1.pointer.child) != .array;
    // if (info1 != .pointer or info1.pointer.size != .slice) {
    //     @compileError("arry error type");
    // }
}

pub fn matrix_init(m: anytype, v: @TypeOf(m[0][0])) void {
    comptime {
        const type1 = type_of_matrix(@TypeOf(m));
        if (type1 == 0 or @typeInfo(@TypeOf(m)) != .pointer or
            @typeInfo(@TypeOf(m)).pointer.is_const)
        {
            @compileError("type is not nultable matrix");
        }
        // const T=@TypeOf(v);
        // if (@TypeOf(m[0][0]) != @TypeOf(v)) {
        //     @compileError("type is not same");
        // }
    }
    for (0..m.len) |i| {
        for (0..m[0].len) |j| {
            m[i][j] = v;
        }
    }
}

/// plus
pub fn matrix_plus(m1: anytype, m2: anytype, re: anytype) void {
    comptime {
        const type1 = type_of_matrix(m1);
        const type2 = type_of_matrix(m2);
        if (type1 == 0 or type2 == 0) {
            @compileError("matrix error type");
        }

        const info = @typeInfo(@TypeOf(re));
        if (info != .pointer or info.pointer.is_const) {
            @compileError("error type");
        }
    }
    if (m1.len != m2.len or m1.len != re.len or m1[0].len != re[0].len or m1[0].len != m2[0].len) {
        return;
    }
    for (0..re.len) |i| {
        for (0..re[0].len) |j| {
            re[i][j] = m1[i][j] + m2[i][j];
        }
    }
}
// // minus
// minus - 矩阵减法
pub fn matrix_minus(m1: anytype, m2: anytype, re: anytype) void {
    comptime {
        const type1 = type_of_matrix(m1);
        const type2 = type_of_matrix(m2);
        if (type1 == 0 or type2 == 0) {
            @compileError("matrix error type");
        }

        const info = @typeInfo(@TypeOf(re));
        if (info != .pointer or info.pointer.is_const) {
            @compileError("error type");
        }
    }
    if (m1.len != m2.len or m1.len != re.len or m1[0].len != re[0].len or m1[0].len != m2[0].len) {
        return;
    }
    for (0..re.len) |i| {
        for (0..re[0].len) |j| {
            re[i][j] = m1[i][j] - m2[i][j];
        }
    }
}

//identity
pub fn matrix_identity(m: anytype) void {
    // if(self.get_cols()!=self.get_rows()){return ;}
    comptime {
        const type1 = type_of_matrix(@TypeOf(m));
        if (type1 == 0) {
            @compileError("matrix error type");
        }
        const info = @typeInfo(@TypeOf(m));
        if (info != .pointer or info.pointer.is_const) {
            @compileError("type is not multable");
        }
    }
    for (0..m.len) |i| {
        for (0..m[0].len) |j| {
            m[i][j] = if (i == j) 1 else 0;
        }
    }
}
// transpose - 矩阵转置
pub fn matrix_transpose(m: anytype, re: anytype) void {
    comptime {
        const type1 = type_of_matrix(@TypeOf(m));
        if (type1 == 0) {
            @compileError("m is not matrix");
        }
        const type2 = type_of_matrix(@TypeOf(re));
        if (type2 == 0 or type2 == 1 or @typeInfo(@TypeOf(re)).pointer.is_const) {
            @compileError("re is not mutable matrix");
        }
    }
    if (m.len != re[0].len or m[0].len != re.len) {
        return;
    }
    for (0..re.len) |i| {
        for (0..re[0].len) |j| {
            re[i][j] = m[j][i];
        }
    }
}
// ///矩阵求逆
pub fn matrix_inverse(m: anytype, re: anytype, allocator: ?std.mem.Allocator) !void {
    const type1 = type_of_matrix(@TypeOf(m));
    if (type1 == 0) {
        @compileError("type is not matrix");
    }
    const T = @TypeOf(m[0][0]);
    if (type1 == 1 or type1 == 2) {
        const rows = comptime m.len;
        const cols = comptime m[0].len;
        // 检查是否为方阵
        if (rows != cols) {
            return lbmath.Math_Compute_Abandon.Math_Compute_Abandon;
        }
        const dim = rows;

        // 创建工作矩阵和单位矩阵
        var data: [dim][dim]T = undefined;
        // var data_b: [dim][dim]T = undefined;
        // 初始化数据矩阵和单位矩阵
        for (&data, 0..) |*row, i| {
            for (row, 0..) |*element, j| {
                element.* = m[i][j];
                re[i][j] = 0;
            }
            re[i][i] = 1;
        }
        // 高斯-约旦消元法（按列选主元）
        for (0..dim) |i| {
            var pivot_col: usize = i;
            var max_val = @abs(data[i][pivot_col]);
            for (i + 1..dim) |j| {
                const v = @abs(data[i][j]);
                if (v > max_val) {
                    max_val = v;
                    pivot_col = j;
                }
            }

            // 交换列
            if (pivot_col != i) {
                // 交换 data 矩阵的列
                for (0..dim) |k| {
                    const temp = data[k][i];
                    data[k][i] = data[k][pivot_col];
                    data[k][pivot_col] = temp;
                    const temp1 = re[k][i];
                    re[k][i] = re[k][pivot_col];
                    re[k][pivot_col] = temp1;
                }
            }

            // 检查主元是否为零
            const pivot = data[i][i];
            // if (math.abs(pivot) <= 1e-11)
            if (lbmath.approxEqAbs(T, pivot, 0, null)) {
                return error.Math_Compute_Abandon;
            }

            // 归一化主列
            for (0..dim) |k| {
                data[k][i] = data[k][i] / pivot;
                re[k][i] = re[k][i] / pivot;
            }

            // 消元
            for (0..dim) |j| {
                if (j == i) continue;
                const factor = data[i][j];
                // 下面是否可以优化,data可以从i到dim,data_b是从0到dim
                for (0..dim) |k| {
                    data[k][j] -= factor * data[k][i];
                    re[k][j] -= factor * re[k][i];
                }
            }
        }
    } else {
        const rows = m.len;
        const cols = m[0].len;
        // 检查是否为方阵
        if (rows != cols) {
            return lbmath.Math_Compute_Abandon.Math_Compute_Abandon;
        }
        const dim = rows;

        // 创建工作矩阵和单位矩阵
        const allocator1 = allocator orelse return lbmath.Math_Compute_Abandon.Math_Compute_Abandon;
        // var data: [dim][dim]T = undefined;
        const data = try allocator1.alloc([]T, dim);
        defer allocator1.free(data);
        for (0..dim) |i| {
            data[i] = try allocator1.alloc(T, dim);
        }
        defer {
            for (0..dim) |i| {
                allocator1.free(data[i]);
            }
        }
        // var data_b: [dim][dim]T = undefined;
        // 初始化数据矩阵和单位矩阵
        for (&data, 0..) |*row, i| {
            for (row, 0..) |*element, j| {
                element.* = m[i][j];
                re[i][j] = 0;
            }
            re[i][i] = 1;
        }
        // 高斯-约旦消元法（按列选主元）
        for (0..dim) |i| {
            var pivot_col: usize = i;
            var max_val = @abs(data[i][pivot_col]);
            for (i + 1..dim) |j| {
                const v = @abs(data[i][j]);
                if (v > max_val) {
                    max_val = v;
                    pivot_col = j;
                }
            }

            // 交换列
            if (pivot_col != i) {
                // 交换 data 矩阵的列
                for (0..dim) |k| {
                    const temp = data[k][i];
                    data[k][i] = data[k][pivot_col];
                    data[k][pivot_col] = temp;
                    const temp1 = re[k][i];
                    re[k][i] = re[k][pivot_col];
                    re[k][pivot_col] = temp1;
                }
            }

            // 检查主元是否为零
            const pivot = data[i][i];
            // if (math.abs(pivot) <= 1e-11)
            if (lbmath.approxEqAbs(T, pivot, 0, null)) {
                return error.Math_Compute_Abandon;
            }

            // 归一化主列
            for (0..dim) |k| {
                data[k][i] = data[k][i] / pivot;
                re[k][i] = re[k][i] / pivot;
            }

            // 消元
            for (0..dim) |j| {
                if (j == i) continue;
                const factor = data[i][j];
                // 下面是否可以优化,data可以从i到dim,data_b是从0到dim
                for (0..dim) |k| {
                    data[k][j] -= factor * data[k][i];
                    re[k][j] -= factor * re[k][i];
                }
            }
        }
    }
}
test matrix_detn {
    var mat: [4][4]f64 = undefined;
    matrix_init(&mat, 1);

    // var mat = Matrix_f64_3x3.init(1);
    // mat.identity();
    mat[0][0] = 2;
    mat[0][1] = 3;
    mat[0][2] = 1;
    mat[0][3] = 5;
    mat[1][0] = 4;
    mat[1][1] = -2;
    mat[1][2] = 0;
    mat[1][3] = 1;
    mat[2][0] = 1;
    mat[2][1] = 0;
    mat[2][2] = 3;
    mat[2][3] = 2;
    mat[3][0] = 0;
    mat[3][1] = 2;
    mat[3][2] = -1;
    mat[3][3] = 4;

    const det = matrix_detn(mat, null);
    const det1 = matrix_detn(&mat, null);

    try std.testing.expectEqual(det, -109);
    try std.testing.expectEqual(det1, -109);

    dprint("det:{}\n", .{det});
    var mati: [4][4]f64 = undefined;

    try matrix_inverse(&mat, &mati, null);
    dprint("mati :{any}\n", .{mati});
    // 0.311926605504587 0.146788990825688 -0.211009174311927 -0.321100917431193
    // 0.504587155963303 -0.174311926605505 -0.311926605504587 -0.431192660550459
    // 0.055045871559999999633 -0.0917431192660551 0.256880733944954 -0.174311926605505
    // -0.238532110091743 0.0642201834862385 0.220183486238532 0.422018348623853

}
// , allocator: ?std.mem.Allocator
pub fn matrix_detn(m: anytype, allocator: ?std.mem.Allocator) @TypeOf(m[0][0]) {
    const type1 = type_of_matrix(@TypeOf(m));
    if (type1 == 0) {
        @compileError("type is not matrix");
    }
    const T = @TypeOf(m[0][0]);
    if (type1 == 1 or type1 == 2) {
        const rows = comptime m.len;
        const cols = comptime m[0].len;
        if (rows != cols) {
            return 0;
        }
        var datan: [rows][cols]T = undefined;
        for (0..rows) |i| {
            for (0..cols) |j| {
                datan[i][j] = m[i][j];
            }
        }

        var swap_num: usize = 0;
        for (0..(rows - 1)) |i| {
            //选择主行
            var temp_k: usize = i;
            for ((i + 1)..rows) |k| {
                if (@abs(datan[k][i]) > @abs(datan[temp_k][i])) {
                    temp_k = k;
                }
            }
            //计算符号变换
            swap_num += (temp_k - i);
            // 交换行
            for (0..cols) |j| {
                const temp = datan[i][j];
                datan[i][j] = datan[temp_k][j];
                datan[temp_k][j] = temp;
            }
            if (lbmath.approxEqAbs(T, datan[i][i], 0, null)) {
                return 0;
            }
            // 行消元
            for ((i + 1)..rows) |k| {
                const fac = datan[k][i] / datan[i][i];
                // 下面是否可以优化i到cols
                for (0..cols) |j| {
                    datan[k][j] -= fac * datan[i][j];
                }
            }
        }
        var re: T = 1.0;
        for (0..rows) |i| {
            re *= datan[i][i];
        }
        const tempt: T = @floatFromInt(swap_num % 2);
        // 符号变换
        re *= (1 - 2 * tempt);
        return re;
    } else {
        const rows = m.len;
        const cols = m[0].len;
        if (rows != cols) {
            return 0;
        }
        const allocator1 = allocator orelse return 0;

        // var datan: [rows][cols]T = undefined;
        const datan = allocator1.alloc([]T, rows) catch return 0;
        defer allocator1.free(datan);
        for (0..rows) |i| {
            datan[i] = allocator1.alloc(T, cols) catch return 0;
        }
        defer {
            for (0..rows) |i| {
                allocator1.free(datan[i]);
            }
        }

        for (0..rows) |i| {
            for (0..cols) |j| {
                datan[i][j] = m[i][j];
            }
        }

        var swap_num: usize = 0;
        for (0..(rows - 1)) |i| {
            //选择主行
            var temp_k: usize = i;
            for ((i + 1)..rows) |k| {
                if (@abs(datan[k][i]) > @abs(datan[temp_k][i])) {
                    temp_k = k;
                }
            }
            //计算符号变换
            swap_num += (temp_k - i);
            // 交换行
            for (0..cols) |j| {
                const temp = datan[i][j];
                datan[i][j] = datan[temp_k][j];
                datan[temp_k][j] = temp;
            }
            if (lbmath.approxEqAbs(T, datan[i][i], 0, null)) {
                return 0;
            }
            // 行消元
            for ((i + 1)..rows) |k| {
                const fac = datan[k][i] / datan[i][i];
                // 下面是否可以优化i到cols
                for (0..cols) |j| {
                    datan[k][j] -= fac * datan[i][j];
                }
            }
        }
        var re: T = 1.0;
        for (0..rows) |i| {
            re *= datan[i][i];
        }
        const tempt: T = @floatFromInt(swap_num % 2);
        // 符号变换
        re *= (1 - 2 * tempt);
        return re;
    }
    return 0;
}

test scala_mult_matrix {
    var m: [3][2]f32 = undefined;
    var re: [3][2]f32 = undefined;
    // const T = @TypeOf(&m);
    // if (type_is_vector(T)) {
    //     @compileLog("m is vector*****##");
    //     if (@typeInfo(@typeInfo(@typeInfo(@TypeOf(&m)).pointer.child).array.child) != .array) {
    //         @compileLog("is array ***");
    //     }
    // }
    matrix_init(&m, 1);

    scala_mult_matrix(2.0, m, &re);
    std.debug.print("scala mult matrix:{}\n", .{re[0][0]});

    const m1 = try std.testing.allocator.alloc([]f32, 3);
    defer std.testing.allocator.free(m1);
    for (0..m1.len) |i| {
        m1[i] = try std.testing.allocator.alloc(f32, 2);
    }
    defer {
        for (0..m1.len) |i| {
            std.testing.allocator.free(m1[i]);
        }
    }
    // matrix_init(m1, 1);
    matrix_identity(m1);
    const re1 = try std.testing.allocator.alloc([]f32, 3);
    defer std.testing.allocator.free(re1);
    for (0..re1.len) |i| {
        re1[i] = try std.testing.allocator.alloc(f32, 2);
    }
    defer {
        for (0..re1.len) |i| {
            std.testing.allocator.free(re1[i]);
        }
    }
    scala_mult_matrix(-2.1, m1, re1);
    std.debug.print("scala mult matrix:{}\n", .{re1[0][0]});
    var arr1: [4]i32 = .{ -1, 2, 3, 3 };
    scala_mult_matrix(1, arr1, &arr1);
    std.debug.print("scala mult vector:{}\n", .{arr1[0]});
}
pub fn scala_mult_matrix(fac: anytype, m: anytype, re: anytype) void {
    // const type1= type_of_matrix(@TypeOf(m));
    if (@typeInfo(@TypeOf(re)) != .pointer or
        @typeInfo(@TypeOf(re)).pointer.is_const)
    {
        @compileError("re is not multable");
    }
    // const info = @typeInfo(@TypeOf(m));
    const type1 = type_of_matrix(@TypeOf(m));
    if (type_is_vector(@TypeOf(m)) and type_is_vector(@TypeOf(re))) {
        for (0..m.len) |i| {
            re[i] = fac * m[i];
        }
    } else if (type1 == 2 or type1 == 3) {
        for (0..m.len) |i| {
            for (0..m[0].len) |j| {
                re[i][j] = fac * m[i][j];
            }
        }
    }
}
test mult_matrix {
    // const Matrix_f64_4x3 = LBMatrix(f64, 4, 3);
    // const Matrix_f64_3x5 = LBMatrix(f64, 3, 5);
    var m1: [4][3]f64 = undefined;
    var m2: [3][5]f64 = undefined;

    matrix_init(&m1, 1);
    matrix_init(&m2, 2);
    var mat: [4][5]f64 = undefined;
    try mult_matrix(m1, m2, &mat);
    dprint("mat ele:{} {} {}\n", .{ mat[0][0], mat[1][1], mat[2][2] });

    // const mat = mat1.mult(mat2);
    // dprint("lb matrix rows:{} cols:{}\n", .{ @TypeOf(mat).get_rows(), @TypeOf(mat).get_cols() });
    // for (2..1) |i| {
    //     dprint("****************{}\n", .{i});
    // }
    // dprint("lb matrix rows:{} cols:{}\n",.{mat1.get_rows(),mat1.get_cols()});

}

pub fn mult_matrix(
    m1: anytype,
    m2: anytype,
    re: anytype,
) !void {
    comptime {
        const type1 = type_of_matrix(@TypeOf(m1));
        const type2 = type_of_matrix(@TypeOf(m2));
        if (type1 == 0 or type2 == 0) {
            @compileError("error type");
        }
        const info = @typeInfo(@TypeOf(re));
        // 检查是否为可变切片 []T或*[N]T
        const is_mutable_slice =
            info == .pointer and
            (info.pointer.size == .slice or (info.pointer.size == .one or @typeInfo(info.pointer.child) == .array)) and
            !info.pointer.is_const;
        if (!is_mutable_slice) {
            @compileError("is not a mutable slice");
        }
    }

    if (m1.len == 0 or m2.len == 0 or m1[0].len != m2.len) {
        return lbmath.Math_Compute_Abandon.Math_Compute_Abandon;
    }
    const T = @TypeOf(m1[0][0]);

    for (0..m1.len) |i| {
        for (0..m2[0].len) |j| {
            var sum: T = 0;
            for (0..m1[0].len) |k| {
                sum += m1[i][k] * m2[k][j];
            }
            re[i][j] = sum;
        }
    }
}

test matrix_mult_array {
    const m1 = try std.testing.allocator.alloc([]f32, 5);
    defer std.testing.allocator.free(m1);
    for (0..m1.len) |i| {
        m1[i] = try std.testing.allocator.alloc(f32, 3);
    }
    defer {
        for (0..m1.len) |i| {
            std.testing.allocator.free(m1[i]);
        }
    }
    matrix_init(m1, 1);
    const arr1 = try std.testing.allocator.alloc(f32, 3);
    defer std.testing.allocator.free(arr1);
    @memset(arr1, 0.5);
    // for (0..arr1.len) |i| {
    //     arr1[i] = 0.5;
    // }
    var re1: [5]f32 = undefined;
    try matrix_mult_array(m1, arr1, &re1);
    std.debug.print("matrix mult array{} {} {}\n", .{ re1[0], re1[1], re1[2] });
    var m2: [5][2]f32 = .{ .{ 0, 0 }, .{ 0, 0 }, .{ 0, 0 }, .{ 0, 0 }, .{ 0, 0 } };
    matrix_init(&m2, 0.1);
    const arr2: [2]f32 = .{ -0.4, 2 };
    const re2 = try std.testing.allocator.alloc(f32, 5);
    defer std.testing.allocator.free(re2);
    try matrix_mult_array(m2, arr2, re2);
    std.debug.print("matrix mult array{} {} {}\n", .{ re2[0], re2[1], re2[2] });
}

pub fn matrix_mult_array(
    m1: anytype,
    arr: anytype,
    re: anytype,
) !void {
    comptime {
        const type1 = type_of_matrix(@TypeOf(m1));
        if (type1 == 0) {
            @compileError("error type");
        }
        // const info1 = @typeInfo(@TypeOf(arr));
        // if (info1 != .pointer or info1.pointer.size != .slice) {
        //     @compileError("arry error type");
        // }
        if (!type_is_vector(@TypeOf(arr))) {
            @compileError("arry error type");
        }
        const info = @typeInfo(@TypeOf(re));
        // 检查是否为可变切片 []T或*[N]T
        const is_mutable_slice =
            info == .pointer and
            (info.pointer.size == .slice or (info.pointer.size == .one or
                @typeInfo(info.pointer.child) == .array)) and
            !info.pointer.is_const;
        if (!is_mutable_slice) {
            @compileError("is not a mutable slice");
        }
    }
    if (m1.len == 0 or m1[0].len != arr.len) {
        return lbmath.Math_Compute_Abandon.Math_Compute_Abandon;
    }
    const T = @TypeOf(m1[0][0]);

    for (0..m1.len) |i| {
        var sum: T = 0;
        for (0..arr.len) |j| {
            sum += m1[i][j] * arr[j];
        }
        re[i] = sum;
    }
}
pub fn array_mult_matrix(
    arr: anytype,
    m1: anytype,
    re: anytype,
) !void {
    comptime {
        const type1 = type_of_matrix(@TypeOf(m1));
        if (type1 == 0) {
            @compileError("error type");
        }
        // const info1 = @typeInfo(@TypeOf(arr));
        // if (info1 != .pointer or info1.pointer.size != .slice) {
        //     @compileError("arry error type");
        // }
        if (!type_is_vector(@TypeOf(arr))) {
            @compileError("arry error type");
        }
        const info = @typeInfo(@TypeOf(re));
        // 检查是否为可变切片 []T或*[N]T
        const is_mutable_slice =
            info == .pointer and
            (info.pointer.size == .slice or (info.pointer.size == .one or
                @typeInfo(info.pointer.child) == .array)) and
            !info.pointer.is_const;
        if (!is_mutable_slice) {
            @compileError("is not a mutable slice");
        }
    }
    if (arr.len != m1.len or m1.len == 0) {
        return lbmath.Math_Compute_Abandon.Math_Compute_Abandon;
    }
    const T = @TypeOf(m1[0][0]);
    for (0..m1[0].len) |i| {
        var sum: T = 0;
        for (0..arr.len) |j| {
            sum += m1[j][i] * arr[j];
        }
        re[i] = sum;
    }
}

// 原地 LUP 分解（部分主元），覆盖 a 为 LU，perm 长度为 n，返回交换次数
pub fn lupDecompose(a: anytype, perm: []usize) !usize {
    comptime {
        const type1 = type_of_matrix(@TypeOf(a));
        if (type1 != 2 and type1 != 3) {
            @compileError("is not a matrix");
        }
        if (@typeInfo(@TypeOf(a)).pointer.is_const) {
            @compileError("is not a mutable matrix");
        }
    }
    const T = @TypeOf(a[0][0]);

    var swaps: usize = 0;
    const n = perm.len;
    // 初始化 perm（perm[i] 表示原矩阵的哪一行现在位于第 i 行）
    for (0..n) |i| {
        perm[i] = i;
    }

    for (0..n) |k| {
        // 找当前列 k 的最大主元（绝对值）
        var maxRow = k;
        var maxVal = @abs(a[k][k]);
        for ((k + 1)..n) |i| {
            const v = @abs(a[i][k]);
            if (v > maxVal) {
                maxVal = v;
                maxRow = i;
            }
        }

        if (lbmath.approxEqAbs(T, maxVal, 0, null)) {
            return lbmath.Math_Compute_Abandon.Math_Compute_Abandon;
        }

        if (maxRow != k) {
            // 交换整行 k 与 maxRow
            for (0..n) |j| {
                const t = a[k][j];
                a[k][j] = a[maxRow][j];
                a[maxRow][j] = t;
            }
            // 交换 perm 的对应项
            const tp = perm[k];
            perm[k] = perm[maxRow];
            perm[maxRow] = tp;
            swaps += 1;
        }

        // 消元：把乘子存入下三角
        // i = k + 1;
        const akk = a[k][k];
        // while (i < n) : (i += 1)
        for ((k + 1)..n) |i| {
            const aik = a[i][k] / akk;
            a[i][k] = aik; // 存 L(i,k)
            var j2 = k + 1;
            while (j2 < n) : (j2 += 1) {
                a[i][j2] -= aik * a[k][j2];
            }
        }
    }

    return swaps;
}

/// 对角化增广矩阵
/// 可以用来计算线程方程组ax=b
/// inverse_a可以用来储存变换矩阵
/// x 解空间,为true是解变量
/// 返回是否存在"0=非零"的矛盾，即是否有解
/// 如果b时null，则表示0向量
fn diagonalizable_matrix(
    comptime T: type,
    a: [][]T,
    b: ?[]T,
    inverse_a: ?[][]T,
    x: ?[]bool,
    // allocator: std.mem.Allocator,
) bool {
    const n = a.len;
    // var re = try allocator.alloc([]f64, n);
    if (b) |bv| {
        std.debug.assert(bv.len == n);
        // if (bv.len != n) {
        //     return;
        // }
    }

    // 初始化数据矩阵和单位矩阵
    if (inverse_a) |re_v| {
        for (0..n) |i| {
            // re[i] = try allocator.alloc(f64, n);
            @memset(re_v[i], 0);
            re_v[i][i] = 1;
        }
    }
    if (x) |xv| {
        // if (xv.len != a[0].len) {
        //     return;
        // }
        std.debug.assert(xv.len == a[0].len);
        @memset(xv, true);
    }

    const m = a[0].len;
    const minn = if (m < n) m else n;
    var ii = 0;
    for (0..minn) |i| {
        // 选择行主元
        // 找到第i列中绝对值最大的元素
        var max_row = ii;
        var max_val = @abs(a[ii][i]);
        for (ii + 1..n) |k| {
            const val = @abs(a[k][i]);
            if (val > max_val) {
                max_row = k;
                max_val = val;
            }
        }
        if (lbmath.approxEqAbs(f64, max_val, 0, null)) {
            continue;
        }
        if (x) |xv| {
            xv[i] = false;
        }
        // 交换行
        if (max_row != ii) {
            for (i..m) |j| {
                const temp = a[ii][j];
                a[ii][j] = a[max_row][j];
                a[max_row][j] = temp;
            }
            if (inverse_a) |re_v| {
                for (0..n) |j| {
                    const temp1 = re_v[ii][j];
                    re_v[ii][j] = re_v[max_row][j];
                    re_v[max_row][j] = temp1;
                }
            }

            if (b) |bv| {
                const temp3 = bv[ii];
                bv[ii] = bv[max_row];
                bv[max_row] = temp3;
            }
        }
        // 新增：归一化主元行
        const pivot_val = a[ii][i];
        // 从主元位置开始归一化即可，前面的已经是0
        for (i..m) |j| {
            a[ii][j] /= pivot_val;
        }
        if (inverse_a) |re_v| {
            for (0..n) |j| { // re 是 n x n 的
                re_v[ii][j] /= pivot_val;
            }
        }

        if (b) |bv| {
            bv[ii] /= pivot_val;
        }
        //消元
        for (0..n) |k| {
            if (k == ii) continue;
            const factor = a[k][i];
            for (i..m) |j| {
                a[k][j] -= factor * a[ii][j];
            }
            if (inverse_a) |re_v| {
                for (0..n) |j| {
                    re_v[k][j] -= factor * re_v[ii][j];
                }
            }

            if (b) |bv| {
                bv[k] -= factor * bv[ii];
            }
        }
        ii += 1;
    }

    if (b) |bv| {
        while (ii < n) {
            if (!lbmath.approxEqAbs(T, bv[ii], 0, null)) {
                // if(bv[ii]!=0){
                return false;
            }
            ii += 1;
        }
    }
    return true;
    // return re;
}

/// 根据上面的对角化求解变量
/// x=re*x_represent+re_b
fn solve_Ax_b(
    comptime T: type,
    a: [][]T,
    b: ?[]T,
    x: []bool,
    re: [][]T,
    re_b: []T,
) void {
    // var num: usize = 0;
    // for (x) |xv| {
    //     if (xv) {
    //         num += 0;
    //     }
    // }
    for (0..x.len) |j| {
        re_b[j] = 0;
        for (re[j]) |*v| {
            v.* = 0;
        }
    }
    var i: usize = 0;
    for (0..x.len) |j| {
        if (x[j]) {
            re[j][i] = 1;
            i += 1;
        } else {
            re_b[j] = if (b) |b_vec| b_vec[j - i] else 0;
            // j - i;
            var k: usize = 0;
            for ((j + 1)..x.len) |l| {
                if (x[l]) {
                    re[j][i + k] = -a[j - i][l];
                    k += 1;
                }
            }
        }
    }
}
